## DiD analysis for presidential electoral results 
## Requires a single data file: dataset-mcmv-by-preselectoral-period.RData, 
## Data file was assembled in Assemble_ElectoralPeriod_Pres.R
## Assumes data file is in subfolder "Data"
## Saves tables in .tex format to subfolder "Table"
## Save figures in .pdf to subfolder "Figures"

## Code consists of the following parts:
## 1) minor recoding of variables for analysis
## 2) analysis by cycle (three cycles)
## 3) regressions on combined (stacked) data from all cycles
## 4) assemble results (table and descriptives) reported in main body of the paper
## 5) assemble tables and descriptives for appendices

rm(list=ls())
library(estimatr)
library(plm)
require(lmtest)
library(texreg)
library(MatchIt)
library(xtable)
#library(remotes)
#remotes::install_github("pedrohcgs/DRDID")
library(DRDID)

load("Data/dataset-mcmv-by-preselectoral-period.RData")
cat(comment(d))

## (1) Recoding some variables #####
d$inc.vs <- d$vs  # Rename incumbent party vote share variable 
d$N.pc <- d$N/(d$pop/100) #benefic. per 100 people
d$N.hd <- d$N/(d$hdtot.2010) #beneficiaries relative to housing deficit
d$N.bin <- d$N>0 # Create dichotomous MCMV indicator (main treatment variable)
table(Period=d$elec.period,d$N.bin)

## (2) Two period DiD by cycles #####

## (2a) 2006-2010 electoral cycle 
d1 <- subset(d,elec.period==0|elec.period==1)
  # Prepare data for "traditional" DiD analysis (equivalent to TWFE)
  d1$MCMV <- is.element(d1$Cod_IBGE,d1$Cod_IBGE[which(d1$N.bin==T)])
  d1$POST <- as.numeric(d1$elec.period==1)
  #For non binary treatments, repeat the "treatment" value in the pre-treatment period
  tmp <- subset(d1,select=c(Cod_IBGE,N.pc,N.hd),elec.period==1)
  d1 <- merge(d1,tmp,by="Cod_IBGE",suffixes=c("",".post"),all.x=T)
  #Actual regressions 
  reg1.1 <- lm_robust(inc.vs~POST*MCMV,data=d1, clusters = Cod_IBGE)
  reg1.2 <- lm_robust(inc.vs~POST*N.pc.post,data=d1, clusters = Cod_IBGE)
  reg1.3 <- lm_robust(inc.vs~POST*N.hd.post,data=d1, clusters = Cod_IBGE)

  #Matched analysis
  to.match <- na.omit(subset(d1,elec.period==0
                     ,select=c(Cod_IBGE,MCMV,pib_pc,pop,hdtotrel.2010,mr.2010,reg)))
  m1 <- matchit(MCMV~pib_pc+pop+hdtotrel.2010+I(hdtotrel.2010^2)+I(log(pib_pc))+mr.2010
                ,replace=T,ratio=2,data=to.match)
  m1 <- subset(get_matches(m1),select=c(Cod_IBGE,id,weights))
  dm1 <- merge(m1,d1,by="Cod_IBGE",all.x=T)
  
  regm1.1 <- lm_robust(inc.vs~POST*MCMV,data=dm1, clusters = Cod_IBGE,weights=weights)
  regm1.2 <- lm_robust(inc.vs~POST*N.pc.post,data=dm1, clusters = Cod_IBGE,weights=weights)
  regm1.3 <- lm_robust(inc.vs~POST*N.hd.post,data=dm1, clusters = Cod_IBGE,weights=weights)

  #Conditioning on observables by improved locally efficient DR DID:
  #Need to eliminate NAs and use baseline covariates only
  d1_did <- d1
  d1_did$MCMVn <- as.numeric(d1_did$MCMV)
  d1_did$Cod_IBGEn <- as.numeric(paste(d1_did$Cod_IBGE,"1",sep=""))
  #Generate gdp_pc and pop at baseline (t0) 
    t0 <- subset(d1_did,POST==0,select=c(Cod_IBGE,pop,pib_pc))
    names(t0)[-1] <- paste(names(t0)[-1],"t0",sep="_")
    d1_did <- merge(d1_did,t0,by="Cod_IBGE",all.x=T)
    rm(t0)
  # Get rid of municipalities with missing covariate data
    to.drop <- unique(d1_did$Cod_IBGE[which(is.na(d1_did$pop_t0)|
                                              is.na(d1_did$pib_pc_t0)|
                                              is.na(d1_did$hdtot.2010)|
                                              is.na(d1_did$mr.2010))])
    d1_did <- d1_did[-which(is.element(d1_did$Cod_IBGE,to.drop)),]
    cat("Keeping",length(unique(d1_did$Cod_IBGE)),"out of",length(unique(d1$Cod_IBGE))
    ,"municipalities in the drdid estimation of cycle 1.\n")
  regdid1.1 <- drdid(yname = "inc.vs"
               , tname = "POST"
               , idname = "Cod_IBGEn"
               , dname = "MCMVn"
               , xformla= ~ pop_t0+pib_pc_t0+hdtotrel.2010+I(hdtotrel.2010^2)+log(pib_pc_t0)+mr.2010
               , data = d1_did, panel = TRUE)

  # Robustness: keep only "not yet treated" (discard never treated)
  d1r <- subset(d1,is.na(first.treated)==F)
  regr1.1 <- lm_robust(inc.vs~POST*MCMV,data=d1r, clusters = Cod_IBGE)
  regr1.2 <- lm_robust(inc.vs~POST*N.pc.post,data=d1r, clusters = Cod_IBGE)
  regr1.3 <- lm_robust(inc.vs~POST*N.hd.post,data=d1r, clusters = Cod_IBGE)

# (2b) 2010-2014 electoral cycle
to.keep <- d$Cod_IBGE[which(d$elec.period==1&d$N.bin==F)]#places without MCMV in ep=1
d2 <- subset(d,is.element(Cod_IBGE,to.keep)&(elec.period==1|elec.period==2))
  # Prepare data for "traditional" DiD analysis (equivalent to TWFE)
  d2$MCMV <- is.element(d2$Cod_IBGE,d2$Cod_IBGE[which(d2$N.bin==T)])
  d2$POST <- as.numeric(d2$elec.period==2)
  table(d2$MCMV,d2$elec.period) #check
  #Repeat the "treatment" value in the pre-treatment period
  tmp <- subset(d2,select=c(Cod_IBGE,N.pc,N.hd),elec.period==2)
  d2 <- merge(d2,tmp,by="Cod_IBGE",suffixes=c("",".post"),all.x=T)
  # Actual regressions
  reg2.1 <- lm_robust(inc.vs~POST*MCMV,data=d2, clusters = Cod_IBGE)
  reg2.2 <- lm_robust(inc.vs~POST*N.pc.post,data=d2, clusters = Cod_IBGE)
  reg2.3 <- lm_robust(inc.vs~POST*N.hd.post,data=d2, clusters = Cod_IBGE)

  #match untreated municipalities at start of cycle 2 (period 1)
  to.match <- na.omit(subset(d2,elec.period==1
                           ,select=c(Cod_IBGE,MCMV,pib_pc,pop,hdtotrel.2010,mr.2010,reg)))
  m2 <- matchit(MCMV~pib_pc+pop+hdtotrel.2010+I(hdtotrel.2010^2)+I(log(pib_pc))+mr.2010
              ,replace=T,ratio=2,data=to.match)
  m2 <- subset(get_matches(m2),select=c(Cod_IBGE,id,weights))
  dm2 <- merge(m2,d2,by="Cod_IBGE",all.x=T)
  regm2.1 <- lm_robust(inc.vs~POST*MCMV,data=dm2, clusters = Cod_IBGE, weights=weights)
  regm2.2 <- lm_robust(inc.vs~POST*N.pc.post,data=dm2, clusters = Cod_IBGE, weights=weights)
  regm2.3 <- lm_robust(inc.vs~POST*N.hd.post,data=dm2, clusters = Cod_IBGE, weights=weights)

  #Robustness: conditioning on observables by improved locally efficient DR DID:
  #Need to eliminate NAs and use baseline covariates only
  d2_did <- subset(d2, elec.period==1|elec.period==2)
  d2_did$MCMVn <- as.numeric(d2_did$MCMV)
  d2_did$Cod_IBGEn <- as.numeric(paste(d2_did$Cod_IBGE,"2",sep=""))
  #Generate gdp_pc and pop at baseline (t0) 
  t0 <- subset(d2_did,POST==0,select=c(Cod_IBGE,pop,pib_pc))
  names(t0)[-1] <- paste(names(t0)[-1],"t0",sep="_")
  d2_did <- merge(d2_did,t0,by="Cod_IBGE",all.x=T)
  rm(t0)
  # Get rid of municipalities with missing covariate data
  to.drop <- unique(d2_did$Cod_IBGE[which(is.na(d2_did$pop_t0)|
                                          is.na(d2_did$pib_pc_t0)|
                                          is.na(d2_did$hdtot.2010)|
                                          is.na(d2_did$mr.2010))])
  d2_did <- d2_did[-which(is.element(d2_did$Cod_IBGE,to.drop)),]
  cat("Keeping",length(unique(d2_did$Cod_IBGE)),"out of",length(unique(d2$Cod_IBGE))
    ,"municipalities in the drdid estimation of cycle 2.\n")

  regdid2.1 <- drdid(yname = "inc.vs"
                   , tname = "POST"
                   , idname = "Cod_IBGEn"
                   , dname = "MCMVn"
                   , xformla= ~ pop_t0+pib_pc_t0+hdtotrel.2010+I(hdtotrel.2010^2)+log(pib_pc_t0)+mr.2010
                   , data = d2_did, panel = TRUE)

#Placebo test for electoral period 2 
#Data includes electoral period 0, 1 and 2 
#Only municipalities not treated in period 1
d2p <- subset(d,is.element(Cod_IBGE,to.keep)
              &(elec.period==0|elec.period==1|elec.period==2))
d2p$MCMV <- is.element(d2p$Cod_IBGE,d2p$Cod_IBGE[which(d2p$N.bin==T)])
d2p$POST <- as.numeric(d2p$elec.period==2)
table(d2p$MCMV,d2p$elec.period) #check
reg2p <- lm_robust(inc.vs~as.factor(elec.period)*MCMV,clusters = Cod_IBGE, data=d2p)

# Robustness for period 2
# Keep only "not yet treated" (discard never treated)
d2r <- subset(d2,is.na(first.treated)==F)
reg2r <- lm_robust(inc.vs~POST*MCMV,data=d2r, clusters = Cod_IBGE)

# (2c) 2014-2018 electoral cycle 
to.keep <- intersect(
            d$Cod_IBGE[which(d$elec.period==2&d$N.bin==F)],
            d$Cod_IBGE[which(d$elec.period==1&d$N.bin==F)])
d3 <- subset(d,is.element(Cod_IBGE,to.keep)&(elec.period==2|elec.period==3))
d3$MCMV <- is.element(d3$Cod_IBGE,d3$Cod_IBGE[which(d3$N.bin==T)])
d3$POST <- as.numeric(d3$elec.period==3)
table(d3$MCMV,d3$elec.period) #check
  #Repeat the "treatment" value in the pre-treatment period
  tmp <- subset(d3,select=c(Cod_IBGE,N.pc,N.hd),elec.period==3)
  d3 <- merge(d3,tmp,by="Cod_IBGE",suffixes=c("",".post"),all.x=T)
  #Actual regressions
  reg3.1 <- lm_robust(inc.vs~POST*MCMV,data=d3, clusters = Cod_IBGE)
  reg3.2 <- lm_robust(inc.vs~POST*N.pc.post,data=d3, clusters = Cod_IBGE)
  reg3.3 <- lm_robust(inc.vs~POST*N.hd.post,data=d3, clusters = Cod_IBGE)

  #match untreated municipalities at start of cycle 3 (period 2)
  to.match <- na.omit(subset(d3,elec.period==2
                           ,select=c(Cod_IBGE,MCMV,pib_pc,pop,hdtotrel.2010,reg,mr.2010)))
  m3 <- matchit(MCMV~pib_pc+pop+hdtotrel.2010+I(hdtotrel.2010^2)+I(log(pib_pc))+mr.2010
              ,replace=T,ratio=2,data=to.match)
  m3 <- subset(get_matches(m3),select=c(Cod_IBGE,id,weights))
  dm3 <- merge(m3,d3,by="Cod_IBGE",all.x=T)
  regm3.1 <- lm_robust(inc.vs~POST*MCMV,data=dm3, clusters = Cod_IBGE, weights=weights)
  regm3.2 <- lm_robust(inc.vs~POST*N.pc.post,data=dm3, clusters = Cod_IBGE, weights=weights)
  regm3.3 <- lm_robust(inc.vs~POST*N.hd.post,data=dm3, clusters = Cod_IBGE, weights=weights)

  #Conditioning on observables by improved locally efficient DR DID:
  #Need to eliminate NAs and use baseline covariates only
  d3_did <- d3
  d3_did$MCMVn <- as.numeric(d3_did$MCMV)
  d3_did$Cod_IBGEn <- as.numeric(paste(d3_did$Cod_IBGE,"3",sep=""))
  #Generate gdp_pc and pop at baseline (t0) 
  t0 <- subset(d3_did,POST==0,select=c(Cod_IBGE,pop,pib_pc))
  names(t0)[-1] <- paste(names(t0)[-1],"t0",sep="_")
  d3_did <- merge(d3_did,t0,by="Cod_IBGE",all.x=T)
  rm(t0)
  # Get rid of municipalities with missing covariate data
  to.drop <- unique(d3_did$Cod_IBGE[which(is.na(d3_did$pop_t0)|
                                            is.na(d3_did$pib_pc_t0)|
                                            is.na(d3_did$hdtot.2010)|
                                            is.na(d3_did$mr.2010))])
  d3_did <- d3_did[-which(is.element(d3_did$Cod_IBGE,to.drop)),]
  cat("Keeping",length(unique(d3_did$Cod_IBGE)),"out of",
      length(unique(d3$Cod_IBGE))
      ,"municipalities in the drdid estimation of cycle 3.\n")
  regdid3.1 <- drdid(yname = "inc.vs"
                     , tname = "POST"
                     , idname = "Cod_IBGEn"
                     , dname = "MCMVn"
                     , xformla= ~ pop_t0+pib_pc_t0+hdtotrel.2010+I(hdtotrel.2010^2)+log(pib_pc_t0)+mr.2010
                     , data = d3_did, panel = TRUE)

  #Placebo test for electoral period 3 
  #Data includes electoral period 0, 1, 2, and 3
  #Only municipalities not yet treated by period 2
  d3p <- subset(d,is.element(Cod_IBGE,to.keep)
                &(elec.period==0|elec.period==1|elec.period==2|elec.period==3))
  d3p$MCMV <- is.element(d3p$Cod_IBGE,d3p$Cod_IBGE[which(d3p$N.bin==T)])
  d3p$POST <- as.numeric(d3p$elec.period==3)
  table(d3p$MCMV,d3p$elec.period) #check
  
  reg3p <- lm_robust(inc.vs~as.factor(elec.period)*MCMV,clusters = Cod_IBGE, data=d3p)  
  
## (3) Stack the three cycles for "combined" analyses ###
dpool <- rbind(d1,d2,d3)  
regp.1 <- lm_robust(inc.vs~POST*MCMV,data=dpool, clusters = Cod_IBGE,fixed_effects=~elec.period)
regp.2 <- lm_robust(inc.vs~POST*N.pc.post,data=dpool, clusters = Cod_IBGE,fixed_effects=~elec.period)
regp.3 <- lm_robust(inc.vs~POST*N.hd.post,data=dpool, clusters = Cod_IBGE,fixed_effects=~elec.period)

dmpool <- rbind(dm1,dm2,dm3)
# Can't use weights and fixed effects in lm_robust
# But assembling the matched data set via get_matches() renders weights irrelevant
regmp.1 <- lm_robust(inc.vs~POST*MCMV,data=dmpool, clusters = Cod_IBGE,fixed_effects=~elec.period)
regmp.2 <- lm_robust(inc.vs~POST*N.pc.post,data=dmpool, clusters = Cod_IBGE,fixed_effects=~elec.period)
regmp.3 <- lm_robust(inc.vs~POST*N.hd.post,data=dmpool, clusters = Cod_IBGE,fixed_effects=~elec.period)

# (4) Assemble results in main body of the paper

# (4a) Table 2a - main results for presidential cycles
tab.pres <- texreg(list(b2006.2010=reg1.1,b2010.2014=reg2.1,b2014.2018=reg3.1,bpool=regp.1,bmatched=regmp.1)
, include.ci = FALSE
, caption = "Presidential Elections"
, caption.above = T
, custom.coef.map = list("POST:MCMVTRUE"="MCMV$times$Post")
, digits = 3
, stars = c(0.01, 0.05, 0.1)
, label = "tab:did:pres"
, file = "Tables/tab-2a.tex"
)

# (4b) Mean treatment intensity in treated municipalities, by cycle and overall
mean.N.pc.ep <- round(as.numeric(by(d$N.pc[which(d$N.bin>0)]
                              ,d$elec.period[which(d$N.bin>0)],mean,na.rm=T)),2)
mean.N.pc <- mean(d$N.pc[which(d$N.bin>0)],na.rm=T)

# (4c) Share of votes affected by MCMV 
library(tidyverse)
glimpse(d)
table(d$SIGLA_PARTIDO)
summary(d$QTDE_VOTOS/d$votetot)

#Use actual number of votes in treated municipalities, not just "newly treated" places included in sample
d1 <- d %>% filter(elec.period == 1 & N.bin == 1) 
sum(d1$QTDE_VOTOS,na.rm=T)*coef(reg1.1)["POST:MCMVTRUE"]

d2 <- d %>% filter(elec.period == 2 & N.bin == 1) 
sum(d2$QTDE_VOTOS)*coef(reg2.1)["POST:MCMVTRUE"]

d3 <- d %>% filter(elec.period == 3 & N.bin == 1) 
sum(d3$QTDE_VOTOS)*coef(reg3.1)["POST:MCMVTRUE"]


## (5) DiD tables and stats for the appendix

## (5a) Appendix J: DiD plot Figures

# Cicle 2
table(d2p$MCMV,d2p$elec.period)
pdf(file="Figures/fig-J1a.pdf",
    width=6,height=4)
rt <- lm_robust(vs~as.factor(elec.period)-1,data=subset(d2p,MCMV==T))
rc <-  lm_robust(vs~as.factor(elec.period)-1,data=subset(d2p,MCMV==F))
par(mar=c(3,4,1,5.5))
plot(c(1,3),c(.3,.7),type="n",bty="n",xaxt="n",xlab="",ylab="")
axis(side=1,at=1:3,labels=c(2006,2010,2014))
lines(1:3,summary(rc)$coefficients[,"Estimate"],lty=2)
points(1:3,rc$coefficients,pch=21,bg= "white")
lines(1:2,summary(rt)$coefficients[1:2,"Estimate"],lty=2)
lines(2:3,summary(rt)$coefficients[2:3,"Estimate"],lty=1)
points(1:2,rt$coefficients[1:2],pch=21,bg= "white")
points(3,rt$coefficients[3],pch=21,bg= "black")
axis(side=4,at=c(rc$coefficients[3],
                 rt$coefficients[3]),
     labels=c("Not yet treated\n+\nNever treated","Treat 2"),
     tick=F,las=2,xpd=NA,line=-1.3)
legend(x="top",legend=c("Pre-treatment","Post-treatment"),
       pch=21,pt.bg=c("white",1),lty=c(2,1),bty="n",horiz=T,cex=1.2)
dev.off()

# Cicle 2
table(d3p$MCMV,d3p$elec.period)
pdf(file="Figures/fig-J1b.pdf",
    width=6,height=4)
rt <- lm_robust(vs~as.factor(elec.period)-1,data=subset(d3p,MCMV==T))
rc <-  lm_robust(vs~as.factor(elec.period)-1,data=subset(d3p,MCMV==F))
par(mar=c(3,4,1,5.5))
plot(c(1,4),c(.3,.7),type="n",bty="n",xaxt="n",xlab="",ylab="")
axis(side=1,at=1:4,labels=c(2006,2010,2014,2018))
lines(1:4,summary(rc)$coefficients[,"Estimate"],lty=2)
points(1:4,rc$coefficients,pch=21,bg= "white")
lines(1:3,summary(rt)$coefficients[1:3,"Estimate"],lty=2)
lines(3:4,summary(rt)$coefficients[3:4,"Estimate"],lty=1)
points(1:3,rt$coefficients[1:3],pch=21,bg= "white")
points(4,rt$coefficients[4],pch=21,bg= "black")
axis(side=4,at=c(rc$coefficients[4],
                 rt$coefficients[4]),
     labels=c("Never treated","Treat 3"),
     tick=F,las=2,xpd=NA,line=-1.3)
legend(x="top",legend=c("Pre-treatment","Post-treatment"),
       pch=21,pt.bg=c("white",1),lty=c(2,1),bty="n",horiz=T,cex=1.2)
dev.off()


## (5b) Appendix J: Placebo test Table J1

tab.pres.j1 <- texreg(list(b2010.2014=reg2p,b2014.2018=reg3p)
                        , include.ci = FALSE
                        , caption = "Placebo Tests"
                        , caption.above = T
                        , digits = 3
                        , custom.coef.names = c("2006 Election",
                                                "2010 Election",
                                                "2014 Election",
                                                "MCMV",
                                                "2010 Election times MCMV",
                                                "2014 Election times MCMV",
                                                "2018 Election",
                                                "2018 Election times MCMV")
                        , reorder.coef = c(1,2,3,7,4,5,6,8)
                        , stars = c(0.001, 0.01, 0.05, 0.1)
                        , label = "tab:did:bin:placebo2"
                        , file = "Tables/tab-J1.tex"
                        )

## (5c) Appendix K1: Full set of estimates on raw DiD sample (Table K1)

tab.pres.k1 <- texreg(list(b2006.2010=reg1.1,b2010.2014=reg2.1,b2014.2018=reg3.1,bpool=regp.1)
                   , include.ci = FALSE
                   , caption = "Presidential Elections"
                   , caption.above = T
                   , digits = 3
                   , stars = c(0.01, 0.05, 0.1)
                   , label = "tab:did:pres:full"
                   , file = "Tables/tab-K1.tex"
)

## (5d) Appendix K2: Alternative specifications of MCMV (continuous treatment)

tab.pres.k3 <- texreg(list(pc2006.2010=reg1.2,pc2010.2014=reg2.2,pc2014.2018=reg3.2,pcpool=regp.2,
                           hd2006.2010=reg1.3,hd2010.2014=reg2.3,hd2014.2018=reg3.3,pcpool=regp.3)
                      , include.ci = FALSE
                      , custom.coef.map = list("POST:N.pc.post"="POST:MCMV","POST:N.hd.post"="POST:MCMV")
                      , caption = "Presidential Elections"
                      , caption.above = T
                      , digits = 3
                      , stars = c(0.01, 0.05, 0.1)
                      , label = "tab:did:prescont"
                      , file = "Tables/tab-K3a.tex")


# (5e) Interpreting the effects for models with continuous treatments (in text of Appendix K2)

#Interpreting the pooled model with N.pc
quantile(dpool$N.pc[which(dpool$N.bin>0)],probs=c(.05,.25,.5,.75,.95),na.rm=TRUE)*coef(regp.2)["POST:N.pc.post"]
#Interpreting the pooled model with N.hd
quantile(dpool$N.hd[which(dpool$N.bin>0)],probs=c(.05,.25,.5,.75,.95),na.rm=TRUE)*coef(regp.3)["POST:N.hd.post"]


## (5f) Appendix K3: Conditioning on pre-treatment covariates via matching (Tables K4)

tab.pres.k4 <- texreg(list(m2006.2010=regm1.1,m2010.2014=regm2.1,m2014.2018=regm3.1,mpool=regmp.1)
                      , include.ci = FALSE
                      , caption = "Estimates"
                      , caption.above = T
                      , custom.coef.map = list("POST:MCMVTRUE"="POST:MCMV")
                      , digits = 3
                      , stars = c(0.01, 0.05, 0.1)
                      , label = "tab:did:presmatch"
                      , file = "Tables/tab-K4a.tex")

## (5g) Appendix K3: Conditioning on pre-treatment covariates via DRDID (Table  K5)
## RDDID package generates a strange output, so build table by hand
att = c(regdid1.1$ATT,regdid2.1$ATT,regdid3.1$ATT)
se = c(regdid1.1$se,regdid2.1$se,regdid3.1$se)
n = c(nrow(d1_did),nrow(d2_did),nrow(d3_did))
pvals = 1/2*pt(att/se, df=n-2, lower.tail = FALSE)
stars <- ifelse(pvals<0.01,"***",
                ifelse(pvals>=0.01&pvals<0.05,"**",
                       ifelse(pvals>=0.05&pvals<0.1,"*","")))
the.tab <- data.frame(rbind(att,se,n))
names(the.tab) <- c("d2006.2010","d2010.2014","d2014.2018")
tab.pres.k5 <- xtable(the.tab,digits=3)
caption(tab.pres.k5) <- "Doubly-robust DiD Estimates, Conditioning on Covariates"
label(tab.pres.k5) <- "tab:drdid:pres" 
tab.pres.k5[1,] <- paste("$",sprintf("%01.3f", att),"^{",stars,"}","$",sep="")## replace 1st row 
tab.pres.k5[2,] <- paste("$(",sprintf("%01.3f", se),")$",sep="")## replace 2nd row 
print(tab.pres.k5,digits=3,caption.placement= "top",
      file="Tables/tab-K5a.tex")

## (5h) Appendix K4: Keeping only eventually treated controls (Table K6)  

tab.pres.k6 <- texreg(list(b2010.2014=reg2r)
                      , include.ci = FALSE
                      , caption = "DiD Estimates, keeping only eventually treated controls (2006–2010)"
                      , caption.above = T
                      , digits = 3
                      , stars = c(0.01, 0.05, 0.1)
                      , label = "tab:did:bin:robust1"
                      , file = "Tables/tab-K6.tex"
)
